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1 Introduction and Motivation 

In the 19th century the American mathematician G.W. Hill devised a simple and useful 
approximation for the motion of the moon around the earth with perturbations by the 
sun. To most dynamical astronomers "Hill's Problem" still means a model for motions in 
the solar system in which two nearby bodies move in nearly circular orbits about another 
much larger body at a great distance. These lectures have, however, been motivated by 
a problem in stellar dynamics. 

Consider a star in a star cluster which is itself in orbit about a galaxy (Figure H). 
The star, cluster and galaxy take the place of the moon, earth and sun, respectively. The 
potentials of the cluster and galaxy are not those of a point mass, and the galactic orbits 
of the star and cluster may be far from circular. Nevertheless Hill's problem is a good 
starting point, and it can be modified easily to accommodate the differences. In section 2 
we outline a derivation of Hill's equations, and in section 3 we summarise the appropriate 
extensions. 

Stars gradually escape from star clusters. This has been expected on theoretical 
grounds for many years, ever since a paper by Ambartsumian (1938). Recently, deep 
observations have confirmed this (e.g. Leon et al 2000), by revealing faint streams of stars 
around a number of the globular clusters of our Galaxy. 

Loosely speaking we can say that a star can only escape if its energy exceeds some 
critical energy. The energies of stars change slightly as a result of two-body gravitational 
encounters within clusters, though the time scale on which this happens (the relaxation 
time scale) is very long, of order lO^yr. But the orbital motions of stars within clusters 
have much smaller time scales of order lO^yr, and until recently it was thought that 
escaping stars would leave on a similar time scale. With this assumption, relaxation is 
the bottleneck, and so the escape time scale (e.g. the time taken for half the stars to 
escape) should vary with the relaxation time. 

Nowadays it is possible to simulate the evolution of modest-sized star clusters with 
3 X 10^ or more members, and the predicted escape time scale can be checked empirically. 
Unfortunately the results contradict the theory (Figure |T]). As these simulations require 



considerable extrapolation in particle number to be applicable to real clusters (for 
which ~ 10^) the error of the theory is serious. 

It turns out that the assumption of rapid escape is the main source of error (Fukushige 
& Heggie 2000, Baumgardt 2000a,b). In fact some stars above the escape energy never 
escape (unless some other dynamical process comes into play), and others take much 
longer to escape than had been generally thought. 
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Figure 1. Results of numerical experiments (Aarseth & Heggie, unpublished) on the 
escape of stars from star clusters. The time for half the stars to escape is plotted against the 
original number N of stars in the simulation. Points are averages over several simulations 
at each N , except the largest value. The continuous line shows the prediction of theory, 
i.e. proportional to the relaxation time (see text), and the dashed line is an empirical fit. 

With this motivation, the remaining sections of these lectures are devoted to the 
dynamics of escape. Section 4 analyses the very definition of escape, which is not as 
straightforward as in more familiar situations. The last two sections show some ways 
in which the computation of the escape rate can be approached. The main result of 
section 5 concerns the way in which the time scale of escape depends on the energy, and 
outlines how this resolves the problem of Figure |I[ Much more difficult, from a theoretical 
point of view, is determining the distribution of escape times, and some relevant ideas are 
introduced in Section 6. 



2 Equations of Motion 



2.1 Derivation 



We now outline a derivation of the equations of Hill's problem in the stellar dynamics 
context. To simplify matters as much as possible, however, we treat the cluster and galaxy 
as point masses M^. and Mg ^ M^. (Figure |^), and consider motion of a massless star in 
the same plane of motion. 

If X, y are the coordinates of the star in a rotating frame centred at the cluster centre, 
its velocity relative to the galaxy is (x — uy, y + uj[R + x]). Therefore the Lagrangian for 

If/ \2 / / T-, \\2^ GM„ GMc r, r, r, 

its motion IS £ = - |(x — uy) + [y + uj[R + x)) | H — — — I , where r = x + y 

and B!^ = {R + x)"^ + y"^. (Note here that we are neglecting the motion of the galaxy, 
which will not affect the final approximate set of equations of motion.) 
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Figure 2. Derivation of Hill's equations. The cluster is treated as a point mass Mc in 
uniform circular motion of angular velocity u at a distance R from a point-mass galaxy 
Mg. 

For reasons that will become clear later we switch to a canonical formulation. The 
momenta conjugate to x and y are 

= C± = X - uy 
Py = Cy =y + u{R + x), 

and the Hamiltonian is 



H = xpx + ypy - jC 
I 

= 2^pI + pI) + ^iyPx-[R + x]py) 



GM„ GMr 



R' 



The next step is common to apphcations in the solar system and stellar dynamics 
but has a different name. In applications to the earth-moon-sun problem it is referred 
to as "neglect of the parallax", while in stellar dynamics it is always called a "tidal 
approximation". (Even that phrase betrays how much the subject of stellar dynamics 
owes to the celestial mechanics of the earth-moon-sun system!) We suppose r <^ R and 

approximate ^ — (1 — + - ^ ^ )/-^- We drop constant terms, substitute cu^ — 
R' R 2 R^ 

GMg/R? from the equations of circular motion (again assuming Mc <^ Mg), and replace 
Py ^ Py-\- ojR- (If the other variables are not changed this transformation is canonical.) 
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Next we write down Hamilton's equations 



Then we get H = -(pI+pD + u;(yp^ - xpy) - -u;^(2x^ - y^) 



X = Hp,,Px = -Hx, (1) 
etc. Finally, on eliminating p^ and Py, we get 

X — 2u}y — 600 X — —X (2) 

y + 2ujx = —y, (3) 

which differ from Hill's equations only in notation, and then only slightly. 



2.2 A Generalised Leapfrog 

The leapfrog is a favourite integration algorithm for equations of motion in stellar dynam- 
ics. It is identical to the Verlct algorithm of molecular dynamics. For a one-dimensional 
problem with Hamiltonian p^/2 + V{x), for example, it may be written 

Xn+l ^ Xn + hpn (4) 
Pn+1 = Pn-hV'{Xn+l), (5) 

where h is the time step. Note that the new coordinate is used immediately, which is where 
the algorithm differs from an Euler algorithm. The effect is dramatic, as the long-term 
behaviour of the leapfrog is much better. 

One of the nice properties of the leapfrog is that it is symplectic, like a good Hamil- 
tonian flow. Here we show how to construct a similar algorithm for the Hamiltonian of 
Hill's problem. 

Euler's algorithm would be 

X„+i = X„ hHp{Xn, Pn), Pn+1 = Pn " hH:^{Xn, Pn), 

where we have written x = {x,y) and p = ipx,Py)- We can make this symplectic by 
replacing p„ by Pn+i in the arguments of the derivatives of H, because it then takes the 
form 

Xn+l — ^ p(Xn, Pn+l), Vn — ^ x:(Xn, Pn+l), 

where the generating function T — x„.p„+i -|- hHiyin-, Pn+i)- 



Writing out these equations explicitly for the Hamiltonian of Hill's problem, we obtain 
the algorithm 

+ ^Vn) (6) 
yn+1 = yn + h{Py,n+l " UJXn) (7) 

Px,n = Px,n+1 + h{-UJPy^n+l " 2uj'^Xn + (8) 
Py,n = Py,n+1 + h{^Px,n+l + ^'^Vn + ^^Vn)- (9) 

These equations look horribly implicit, a common difficulty with elementary derivations 
of symplectic methods, but in fact eqs. (||) and (|^) are easily solved explicitly for Pn+i and 
then eqs.@ and (|^) give Xn+i. 



2.3 Elementary Properties 

1. The Hamiltonian Ti is time-independent, and so its value is conserved. Rewrit- 
ing the momenta in terms of the velocity components one finds that this value is 

E = ~^ " 2^^*^^ ^^^'"^ often referred to as the "energy". Again 

there is another name in the celestial mechanics community, who refer to the "Ja- 

cobi constant" C = —2E. In stellar dynamics this term is often applied to E. 

At any rate, one implication is that the motion is bounded to the region in which 
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— -cu^x^ < E. The boundaries of these regions are called HilFs curves 

(Figure!). 

It is sometimes tempting to refer to the expression for E as the Hamiltonian, and 
indeed the right-hand side has the same value as Ti. It is, however, impossible to 
recover the equations of motion from the expression for E. 

2. Hill's equations have two equilibrium solutions, at {x,y) = ±(rt,0), where = 
GMc/(3co'^). In stellar dynamics rt is called the tidal radius or Jacobi radius, and 
in all subjects these points are referred to as the Lagrange points Li and L2. 

3. Hill's equations have an obvious symmetry: if {x{t),y{t)) is a solution, then so is 
{x' (t) , y' (t)) = t), —y{—t)). This is quite useful for studying asymptotic orbits. 
For example, if an orbit tends to Li as t — ^> 00, then the orbit obtained by this 
symmetry tends to Li as t ^ —00. Also, if i;(0) = ?/(0) = 0, the two orbits are the 
same, as they satisfy the same initial conditions. This helps to explain the amount 
of attention that has been paid in the literature to such orbits. 



3 Variants of the Problem 

1. It is not necessary that one of the bodies is massless. Hill's equations are also 
applicable to the relative motion of the moon and earth, under solar perturbation, 
as in Hill's original research. A relatively accessible account of this research is 
Plummer (1918). A modern application is binary asteroids (e.g. Chauvineau & 
Mignard 1990). 
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Figure 3. Hill's curves. 



2. It is not necessary tliat tlie two small bodies are bound. Another application is to 
near-conjunctions of coorbitals (e.g. Murray & Dermott 1999). 

3. When the smallest body is treated as massless, as in the star cluster application, 
it is not necessary that the other bodies are treated as point masses. For a spher- 
ically symmetric galaxy potential and an arbitrary cluster potential 0c the three- 
dimensional equations of motion are 



X — 2ujy + {k^ — 4:Lj'^)x 



y + 2ujx 



Z + UJ z 



d(f)c 
dx 



dy 



dz 



where k is the epicyclic frequency (Chandrasekhar 1942, Binney & Tremaine 1987) 
and the plane of motion of the cluster is the x, y plane. For a point-mass galaxy 
K, = UJ and the previous equations are recovered (when 0^ = —GMc/r.) 

Very often the cluster potential 0c would be chosen to be that of a King model 
(cf. Binney & Tremaine 1987). Qualitatively the most important difference from 
the point-mass potential is that the depth of the potential well is finite. Figure |] 
illustrates the differences which a change of potential can make. 



4. Returning to the point mass case, we now consider the situation in which the motion 
of the cluster is elliptic, with eccentricity e. There is now a formulation using 
the same coordinates as in Figure |^ but scaled by R (so-called rotating, pulsating 



coordinates x, y). For coplanar motion the equations are 



X — 2y 




Mc X 
WgP 



y" + 2y' 



1 . M,y 



1 + e cos / Mg r 



where / denotes differentiation with respect to /, the true anomaly of the cluster 
orbit, and = SP' These equations can be easily derived from the correspond- 
ing formulation of the elliptic restricted problem (Szebehely 1967). An important 
difference from the circular case is that the Hamiltonian is no longer autonomous, 
and there is no Jacobi integral. 

5. One can equally well treat the previous problem in rotating, noB-pulsating coordi- 
nates with origin at the centre of the cluster. For coplanar motion, a point-mass 
galaxy and an arbitrary cluster potential, the equations are 



but the corresponding three-dimensional equations can easily be derived for any 
spherical galaxy potential (Oh, Lin & Aarseth 1992). Here, of course, uj is not 
constant in general. 

6. For a still more general galaxy potential (j)g it is simplest to use non- rotating, non- 
pulsating coordinates, i.e. a coordinate frame with origin at the cluster centre but 
with axes parallel to fixed directions in space. Then the equation of motion takes 
the simple vector form r = — V0c — r.VV(/>g. 

Though this may well be the most useful formulation for non-circular cluster motion, 
and certainly when the potential is not even spherical, one can't help feeling that 
something is lost in this simplicity. For example, in the case of a point mass galaxy 
the equation of motion is 



where R is the unit vector from the galaxy to the cluster. Now the corresponding 
Hamiltonian is time-dependent, and it is no longer obvious that any integral exists. 
But the Jacobi integral is still conserved, taking the form 




50c 



dx 



dy ' 




(10) 



E = -r^ - cj.r X r + -cjV^ cj^(r.R)^ + 0c- 



This is an integral of eq . (|ToD , but not quite an obvious one. 



4 Escape Criteria 



4.1 Escapers 



An escaping star eventually travels far from the cluster, and the cluster potential is 
negligible. If the right sides of eqs.(0) and are neglected we have the approximate 
solution 

X = X + acos{t + (j)) (11) 

y = Yo-^Xt-2asm{t + ^), (12) 

where X, Yq, a and (p are constants, and we have scaled t so that u = 1. Typical orbits 
are shown in Figure ^ even in places where the cluster potential would not be negligible. 
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Figure 4. Orbits in Hill's problem when the cluster potential is neglected. Note that the 
axes are orientated unconventionally. 

Notice that stars like to revolve or spiral anticlockwise, while the axes are such (Figure 
H) that the galaxy is far away at the bottom and the direction of motion of the cluster is 
to the right. Thus the stellar motions are retrograde. This is typical of epicycles, as these 
motions are often called. 

Two orbits drawn in Figure ^ are centred at the location of the cluster. There is a 
family of such orbits, for varying a. When the cluster potential is restored this family be- 
comes a family of stable retrograde satellites of the cluster. Its existence has been known 
for a long time (Jackson 1913, Henon 1969, Benest 1971, Markellos 2000). (In the solar 
system context these are sometimes referred to as eccentric retrograde satellites, but the 
reference to the heliocentric eccentricity is not illuminating in the stellar dynamical con- 
text.) This family, referred to as Family / by Henon, ranges from tiny, almost Keplerian 
orbits around the origin to huge orbits far beyond the tidal radius. As a star cluster loses 
mass by the escape of stars, it is conceivable that some stars in retrograde orbits are left 
on such orbits well outside the tidal radius, and it would be interesting to look for these 
in A^-body simulations. 

Now consider the orbit in Figure ^ passing through the origin. Again such orbits of 
arbitrary size exist (Ross et al 1997). Though severely distorted by the cluster potential 



near the origin, they show that stars can escape, recede to arbitrary distance, and then 
return to the cluster again. Thus distance by itself is no guarantee that escape is perma- 
nent. Rigorous escape criteria can be derived, but, to be frank, in practical terms it is 
quite enough to assume that stars that recede to a few times will escape; the fraction 
of such stars that do not escape is tiny. 



4.2 Non-escapers 

It is easy to obtain a rigorous criterion for non-escape, using the simple idea behind 
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Figure H. A particle at rest at the Li 2 points has energy Ec = -, and any star 

2 n 

with E < Ec, and lying within the Hill curve passing through the Lagrangian points, can 
never escape. 

What now ii E > E^^ We already know one set of orbits on which a star can remain 
inside the cluster forever, even with energy above the escape energy: these are the stable 
retrograde satellites (which move outside the tidal radius only for energies considerably 
above Ec). 

Being stable, these orbits are surrounded by a region of phase space with the same 
properties. This is illustrated by the surfaces of section in Figure ^ Closed invariant 
curves predominate on the left side of the diagrams, which corresponds to retrograde 
motions. At the centre of this nested set of curves is a fixed point corresponding to 
the retrograde periodic orbit. Though these diagrams are plotted for energy E = Ec, 
similar sets of invariant curves are obtained at somewhat higher energies in the standard 
Hill problem (Chauvineau & Mignard 1991, Simo & Stuchi 2000). They correspond to 
retrograde motions of stars permanently remaining inside the cluster and with energies 
above the energy of escape. The chaotic scattering of points on the right-hand half of the 
diagram would, however, generally correspond to escaping orbits for E > Ec- 

It is just possible that such orbits have an astrophysical relevance. In two star clusters 
(Gunn & Griffin 1979, Meylan et al 1991) there are stars whose radial velocity alone 
appears to exceed the escape velocity. Perhaps these are indeed stars permanently bound 
within the cluster at energies above the escape energy. 



5 Escape Rate 

5.1 Motion near the Lagrangian points 

Before attempting to determine the rate at which stars escape, we study orbits in Hill's 
problem a little longer. It is clear from the structure of Hill's curves (Figure |^) that, at 
energies just above the energy of escape, an escaper must make its way at relatively low 
speed through one of the gaps in the Hill curves near Li and L2. Therefore it pays to 
study motions near these points, which can be done by linearisation of the equations of 
motion. 

In the vicinity of {x,y) = {rt,0), when u = GMc = 1, eqs.(^ and (^ take the 
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Figure 5. Surface of section in the two-dimensional Hill's problem at the escape energy 
E = Ec- A point is generated on the surface each time an orbit crosses the line y = 
with y > 0. The edges of the diagram are limiting curves corresponding to the condition 
y = 0. Upper diagram: the potential is that of a model star cluster called a King model 
(from Fukushige & Heggie 2000). Lower diagram: point-mass potential. 



approximate form 



e-277-9e = 
fi + 2^ + 37] = 0, 

where x = rt + ^ and y = f]- These have the general solution 

v) \A-V7) [^-V7J \ (4 + v^)sin(z/^ + ^■ 

(13) 

where A, B, C and 9 are arbitrary constants. On this solution the "energy" is = 
Ec+C2(10y7 + 49) + A5(196 - 40^7). Several cases have interesting properties: 
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1. A = B = C = Q: this is the Lagrange point, where E = (= in general). 

2 rt 

2. B = C = 0: this solution approaches Li as t ^ — oo. It is part of the local unstable 
invariant manifold of Li, and E = Ec- 

3. A = C = 0: this solution approaches Li as t — > oo. It is part of the local stable 
invariant manifold of Li, and E = Ec 

4. A = B = 0: the solution is periodic, and E > Ec- Though derived in a linear 
approximation, there is indeed a family of periodic solutions of the full Hill problem, 
parametrised by E (Liapounov's Theorem, cf. Moser 1968). They are named Family 
a and c (one for each Lagrangian point) in Henon (1969). 

5. A = 0: part of the local stable invariant manifold of the Liapounov orbit (Figure 
I)- 

6. B = 0: part of the local unstable invariant manifold of the Liapounov orbit. 



5.2 The flux of escapers 

Stars escaping from the interior of the star cluster have A < and i? < 0, so that ^ — ±oo 
as t ^ ±oo; thus C^(10\/7 + 49) < E - E^. For fixed E > E^, then, this is stars "inside" 
the tube formed by the stable invariant manifold of the Lagrange point (Figure P). It is 
quite easy to estimate the rate at which the phase space occupied by these escapers flows 
out of the cluster. The general theory is given by MacKay (1990), though some trivial 
generalisation is needed because of the Coriolis forces in Hill's problem. 

The rate of flow of phase space (per unit energy Eq) is 

J" = / x6 {E{x, y, px, Py) - Eq) dydpxdpy, (14) 

J x>0,x=rt 

where the 5-function singles out values of the phase-space variables x, y,Px,Py correspond- 
ing to the required energy. This is readily transformed to 

^ = / t5 (e{^, V, i V) - Eo) dr^dtdf, (15) 



1-1 9 3 
in the notation of Section 5.1. In fact E = -^^ + -t)^ — -C,^ + -^rf + and so 



2'l ^ 2' 

= ^iE,-E,). 

This is a two-dimensional result (i.e. for the coplanar problem). In the three- 
dimensional problem it is found that !F oc [Eo — Ec)"^, with an equally simple coefficient. 
In each case, however, the flux of escaping phase space must be doubled, as there two 
Lagrangian points. 

In order to turn the flux into a time scale for escape, it is only necessary to estimate 
the volume of phase space inside the cluster at energy E. In turn this is given by an 
integral of the form V = J^^rt ^{E — EQ)dxdydpxdpy in two dimensions. This does not 
change much with Eq in the vicinity of the critical energy, and there it is easily seen to 
be 27r times the area inside the last closed Hill curve (Figure 

It follows that the time for escape varies as {E — -Ec)~^ in the three-dimensional 
problem, though this concerns the dimensionless case in which u = 1. When dimen- 
sional factors are reinserted it turns out that the result is a time scale proportional to 

E^ 1 

This is a central result of these lectures. 



{E-E,)^u 



5.3 Numerical Methods 

It is not hard to obtain the rate of escape numerically. One possible procedure is the 
following. 

1. Choose some E > E^, 

2. Select initial conditions at energy E inside the sphere r = rt, according to some 
distribution (cf. Fukushige & Heggie 2000); 

3. Determine the escape time te, defined to be the first time when r > rt {pace the 
problem mentioned in Section 4.1); 

4. Repeat 2-3 many times; 

5. Compute P{t), defined to be the fraction of cases with tf. > t. 

The third item in this procedure requires choice of a numerical integration scheme. 
Many are available, but it is worth mentioning here one of the favourites in this subject, 
which is a fourth-order Hermite scheme (cf. Makino & Aarseth 1992). It is a self-starting 
scheme, and we illustrate it for the one-dimensional equation of motion x = a{x). Suppose 
position and velocity are known at the beginning of a time step of length At, and have 
values xq, vq, respectively. From the equation of motion compute the initial acceleration 
and its initial rate of change, i.e. oq and do, respectively. Compute the predicted position 




Figure 6. Orbits in Hill's problem around one of the Lagrangian points, at a fixed energy 
E just above E^., after Fukushige & Heggie (2000). The potential of the cluster is that of 
a King model. Several orbits are shown which approach a periodic orbit asymptotically. 
Other similar diagrams (for the point mass potential of the usual Hill problem) are given 
in Marchal (1990) and Simo & StucM (2000). 



and velocity at the end of the time step by 



Xp = Xo + voAt + ao— - + do 

2 

Vp = vo + aoAt + oo^-- 

Now compute the acceleration and its derivative at the end of the time step, using Xp and 
Vp. If the results are denoted by ai and hi, respectively, the values of x and v at the end 
of the time step are recomputed by 

At, ^ X 

xi = Xo + —{vo + vi) - —{ai- ao) 

At. , At\. . , 

vi = vo + -^[^0 + ai> - ^{ai - ao>- 



Now we return to the numerical problem of determining the escape rate. A typical 
set of results is shown in Figure |^. Curves at larger t correspond to smaller values of 
E — Ec. It can be seen that these have a horizontal asymptote well above the t-axis; in 
other words, there is a substantial fraction of stars for which the escape time is extremely 
long. This is not unexpected, because of the stable retrograde motions shown in Figure 
^. The fraction of such stars decreases as E increases. We also see, as expected from 



Section 5.2, that the escape times decrease as E increases. Indeed, if we redefine P{t) to 
be the fraction of escapers with escape times t^. > t (i.e. we exclude stars which never 
escape), and if we rescale the values of t hj {E — Ec^, the resulting curves lie very nearly 
together, independent of energy (Fukushige & Heggie 2000). 
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Figure 7. Distribution of escape times from a generalised Hill's problem, for various 
values of the energy, after Fukushige & Heggie (2000). 



5.4 Relaxation and Escape 

In this subsection we now show how the escape rate which we have just determined leads 
to a resolution of the problems with the scaling of A^-body simulations, with which these 
lectures were motivated (Figure |l|). The ideas are based on those given in Baumgardt 
(2000b). 

We imagine that stars are present in a cluster with a distribution of energies n{E). 
This distribution evolves as a result of two processes: (i) relaxation, which is a kind of 
diffusion process with a characteristic time scale tr] and (ii) escape, which takes place on 

Er_ 



a time scale of order P ( — — - ) , where P is the orbital period of the cluster round the 

.E — Ec 



galaxy. 

As a very simple model for this problem we shall consider the toy model defined by 
the differential equation 

dn E^ d^n , 1 (E-E^\^ , ^ 



where n{E)dE is the number of stars with energies in the range {E,E + dE), and the 
Heaviside (unit step) function H confines escape to energies above Ec- 

There are several details missing from this problem. First, in addition to the diffusive 
term (i.e. the first term on the right side) one needs a "drift" term corresponding to 
dynamical friction (cf. Spitzer 1987, for this and other issues in what follows). We have 
also neglected the fact that the coefficient of the diffusion term depends on E and n{E) 
in a complicated way. Next, one needs to take into account the effect on the energies of 
the stars of the slowly changing gravitational potential of the cluster. Finally, we need to 
take into account the stars above the escape energy whose escape time scale appears to 
be infinite. If all those factors were included, we would obtain something close to a full 
Fokker-Planck equation for the evolution of the distribution function in the presence of 
energy-dependent escape. We shall see, however, that this toy model is quite illuminating. 

Let us now scale t by U and let x = {E — Ec)/\Ec\. Then the escape boundary occurs 
at X = 0, and the equation transforms to 

m = d^- """^^^^ ('^^ 

where a = tr/P. Now in astrophysical applications P varies with the crossing time scale 
in a star cluster, and so a varies nearly as N/ log A^, where is the total number of stars 
(cf. Spitzer 1987). Therefore a can be used nearly as a proxy for A^. 

In order to estimate an escape rate we adopt the strategy pioneered by Chandrasekhar 
in this context (Chandrasekhar 1943), which is to look for a separable solution n{x,t) = 
exp{—Xt)y{x), where we expect A > 0. If we impose a no-fiux boundary condition at 
X = —1 (say) and the condition that y{x) ^0 as x oo, then we find that 



Acos |vA(x + 1)| if X < 

B /"^e^^V2|^ _ ,|A/4v/^ - 3/4|^ ^ ,\-X/AV^ - 3/4^, jf ^ > q. 



where A, B are constants. 

While the first of these solutions is elementary, the second deserves some explanation. 
As Maple shows, the solution of the differential equation for y{x) when x > can be 
written in terms of Whittaker functions, and a search through Abramowitz & Stegun 
(1965) shows that these can be expressed as integrals. It is easier, however, to proceed 
directly, though the appropriate methods are not in common use (cf. Burkill 1962). In 
this particular case, for a reason that will become clear, we first change the independent 
variable to z = x^/2. Then the differential equation becomes 

2zy + y + {X~2az)y = 0, (19) 

where a dot denotes a 2;- derivative. Now we get down to business. Motivated by 
the inversion integral for Laplace transforms, we seek a solution in the form y{z) = 

/ exp{sz)f{s)ds, where both the function / and the contour C are to be chosen. Sub- 
Jc 

stituting into eq. ([19|) , we find that we require 



exp{sz)f{s){2zs^ + S + X- 2az)ds = 0. 



No non-trivial choice of / will make the integrand vanish. We can, however, integrate by 
parts to remove the 2;-dependent part of the last factor of the integrand. It follows that 
we require 



[2exp(sz)/(s)(s2 - a)] + J^exp{sz) lf{s){s + A) - ^ {2f{s){s^ - a)}] ds = 



where the first term is the end-point contribution. Now the integral can be made to 
vanish by making the integrand vanish, which in turn requires the solution of a separable 
first-order differential equation. (Without the precaution of changing from x to z, this 
would have been a second-order equation.) This gives the integrand in eq.(|l8|b). To make 
the end-point contribution vanish, we note that we require a function y{z) vanishing as 
z ^ oo, and the exponential factor in the integrand has this behaviour if we restrict the 
contour to Res < 0. One obvious choice for end-point is s = —oo. For the other we 
choose the negative root of /(s)(s^ — a), i.e. s = —y/a, which works if A < 4:^/a. In 
fact the more stringent condition is the integrability of f{s) at this point, which requires 
A < ^/a. 

Now we must match y and y' a.t x = 0. Evaluating the integral at x = is straight- 
forward, and the transformation s = —^/a{l + 2t) gives a standard integral for a beta 
function. In order to evaluate the derivative y'{0+) one cannot simply differentiate the 
integral and substitute x = 0. For one thing the resulting integral diverges as s — >■ —oo. 
This behaviour is killed by the exponential if x is small and positive, and in this case one 
can approximate the other factors in the integrand by their asymptotic form as s ^ — oo. 
Again one obtains a standard integral, this time for a gamma function. 

In the end one finds that the relation between A and a is 

„ / A 3\ 



tanVA = W^ ) ^ f. (20) 



A^/a 4 / 

As A/ ^/a ^1 — , the gamma function in the denominator tends to infinity, and so \/X 0. 
Thus there is an asymptotic regime such that a ^ and A ~ ^/a. If, on the other hand, 
\/y/a — s> 0+, it is clear that the right hand side of eq. (|20|) tends to infinity, and so 
A 7r^/4. Numerical study shows that there is a single solution which joins these two 
asymptotic regimes. 

In the second asymptotic regime {a — > oo, i.e. large A^), escape is very efficient, and 
the time scale for loss of mass, 1/A, is determined by relaxation. Recalling that we have 
scaled time by the relaxation time, it follows that the time to lose half the mass, say, varies 
as tr- In case a is small, however, escape is the bottleneck, and the escape timescale, in 
units of the relaxation time, increases as A^ (or a) decreases. In fact in this regime we 
expect the half mass time to vary nearly as t^/ \/N. Since tr itself varies nearly as A^ (in 
the units of Figure P, we expect a time scale varying as t^^. 

These results correspond qualitatively to what is observed (Figure 0). It should be 
stressed, however, that the value of this toy model is purely qualitative. When one studies 
simulations of the evolution of star clusters quantitatively (Baumgardt 2000b, or those in 
Figure |I|) it is found that, in the case of small A^, the actual scaling is more like tf^'^. 



We now outline Baumgardt's argument which leads to this scaling. We assume that 
the distribution of escapers (i.e. those with E > is nearly in equilibrium. Then 
eq.(|T^ shows that the width of the distribution is x ~ (This scaling can also be 

seen in eq.(|18|b).) The number of such escapers is proportional to this width, and can be 
estimated to be of order Na~^/'^. The escape time scale at this energy is of order l/(ax^), 
~ and therefore the rate of escape is of order iVa^/^. Thus the time scale for losing 

(say) half the mass is of order a~^/* in units of the relaxation time, i.e. the time scale of 
mass loss varies almost as t^/'^. 

It is the assumption that the distribution of escapers reaches equilibrium which dis- 
tinguishes this estimate from the toy model discussed previously, but the reason for this 
difference is not understood. 

6 Distribution of Escape Times 

The results of the previous section relate to the time scale on which stars escape, and 
we conclude with some discussion of the actual distribution of escape times. This issue 
has been studied in a fairly wide variety of problems (e.g. those discussed in the book 
by Wiggins 1992, and Kandrup et al 1999). In some problems the distribution is found 
to be approximately exponential, and in others it is better approximated by a power law. 
For escape in Hill's problem, the numerical experiments summarised in Section 5.3 show 
that the distribution is approximately a power law, over the range of escape times that 
are relevant in applications and amenable to numerical study (Fukushige & Heggie 2000). 

In this section we shall not even come close to obtaining the distribution of escape 
times numerically. We shall, however, introduce two tools which show us how to think 
about this problem. One is a suitable theoretical framework (called turnstile dynamics), 
and the other is a toy model (Henon 1988) which serves two purposes: (i) it can be used 
to illustrate turnstile dynamics, and (ii) it was inspired by Hill's problem. 

6.1 Henon's Toy Model 

We already presented a surface of section for Hill's problem, and Henon's model could have 
been devised with the properties of the corresponding Poincare map in mind. Physically, 
however, it can be thought of as the problem of a ball falling under gravity and bouncing 
off two disks (Figure |^). 

When the radius of the disks is very large, Henon showed that the relation between 
conditions at each bounce takes a particularly simple form, which is 

Xj^i = Xj cosh ip + Wj sinh ip — sj (cosh — 1) 

Wj+i = Xj sinh ip + Wj cosh ip — {sj cosh ip + Sj+i) tanh — , 

where ip is a parameter (related to the radius of disks, the strength of gravity, etc.), Xj is 
the x-coordinate at the jth bounce, Wj is the tangential velocity component at this time, 
and Sj = signXj. (There is a tiny subtlety at Xj = 0, which we ignore in this exposition.) 

The only non-linearities in this problem are the terms with s's. Otherwise the map is 
just a hyperbolic rotation about the point X = ±1, w = 0, in the left and right halves of 
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Figure 8. Henon's billiard model for Hill's problem. 



the X, w plane, as appropriate. It is only when a point moves from one half to the other 
(across the discontinuity in the surface off which the ball bounces) that anything different 
happens. 

These two points are fixed points of the map. As usual in such situations, a fixed point 
corresponds to a periodic motion, which here refers to the ball bouncing repeatedly off 
either of the highest points of the disks (Figure |^). These motions are obviously unstable, 
and the fixed points on the surface of section have local stable and unstable invariant 
manifolds which are segments of the lines X = ±1 ±w (Figure |^). 




Figure 9. Schematic surface of section for Henon's model. The dashed lines are the local 
stable and unstable invariant manifolds of the fixed points at (±1, 0). 

What has this to do with Hill's problem? For one thing the unstable periodic orbits 
have an analogy (in Hill's problem) with the Liapounov orbits mentioned in Sec. 5.1. Using 
the linearised equations derived there it is also possible to derive equations for the local 
stable and unstable invariant manifolds of the corresponding fixed points on the surface 
of section (Figure 0). The main difference between the two models is the absence, in 
Henon's model, of anything comparable with a limiting curve. 
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Figure 10. Outline surface of section for Hill's problem, at some energy E just above 
Ec- The small elliptic arcs are the local stable and unstable manifolds of the fixed points 
corresponding to the Liapounov orbits, and the large curves are the limiting curves. 

6.2 Turnstile Dynamics 

In Figure it is fairly obvious how to define the part of the section "inside" the cluster, 
and how to define the part outside. In Figure ^, despite the absence of limiting curves, 
we shall define the inside and the outside by the naive resemblance of the two pictures. 
To be more precise, the inside (i?i) will be defined as the rhombus lying within the stable 
manifolds of the fixed points, and the outside {R2) as everything else (Figure ^. This at 
least makes clear that the boundary between the two regions is to be defined by pieces 
of the stable and unstable manifolds. This is one of the main procedural points in the 
theory of turnstile dynamics (Wiggins 1992), which we now introduce via this example. 
In order to apply this theory to Hill's problem, there also we would have to define the 
inside and the outside a little more carefully near the fixed points, though we shall not 
dwell on the details here. 

The problem of escape in Henon's model now focuses on the parts of the surface of 
section which, under the Poincare map, are exchanged between regions Ri and R2. A 
direct calculation (Roy 2000) shows that the region which, on one iteration of Henon's 
map, leaves the region Ri consists of the union of two triangles. One of these is shown on 
Figure ^ and labelled Li^2(l)- The notation, which comes from Wiggins (1992), indicates 
that this region is a lobe, which moves from region Ri to region R2 on one iteration. In 
the usual situation considered by Wiggins, a lobe is bounded by parts of the stable and 
unstable manifolds of fixed points. This is only partly true in Henon's model. Two parts 
of the boundary of the little triangle on Figure ^ have this property: the lower right, 
which is part of the unstable manifold of the right-hand fixed point, and the boundary 
at upper right, which consists of part of the stable manifold of the left-hand point. The 
discontinuity at X = provides the remaining part of the boundary. 

We now consider capture of phase space from the region R2 into the region Ri. Again 
we have two triangular regions, one of which is shown in Figure ^, and labelled L2,i(l), 



as the reader should by now appreciate. Also shown in the Figure are successive iterates 
of this lobe under the Henon map H. It can be seen that these remain inside Ri until the 
map is iterated 5 times. The region if'^(L2,i(l)), which is the black triangle furthest to 
the lower right, intersects Li^2(l), and after one further iteration this intersection leaves 
region Ri (and does so forever, actually. The number of iterations that elapse before such 
an intersection takes place depends on the value of ip, of course.) 

Now we can see how the distribution of escape times can be analysed, at least in 
principle. Imagine that, at t = (where t counts the number of iterations) the region Ri is 
filled uniformly with points. At time t = 1, the area occupied by Li_2(l) escapes. The same 
happens at times t = 2, 3, and 4. At time 5, however, the number of points that escape is 
given by the area of Li 2(l)\-f^^(i^2,i(l))- At time 6 the area is now Li^2(l)\(-f^^(-^2,i(l)) U 
/7^(L2,i(l))), and so on. 

Henon's toy problem is unusual in that some of these calculations can be carried out by 
elementary means. In almost all problems, by contrast, the work is necessarily numerical. 
Nevertheless the ideas of turnstile dynamics help to economise the work. The naive way of 
computing a distribution of escape times, as in Sec. 5. 3, is to distribute points throughout 
region Ri and measure how long they take to escape. We now see, however, that we only 
need to consider the dynamics of points within Li.2(l) in order to reach the same results. 
This concentrates the numerical work where it is actually needed. 

When we apply these notions to Hill's problem, a number of additional complicating 
factors arise. In the first place the area on the surface of section is not proportional to the 
volume of phase space (Binney et al 1985), and therefore does not yield an appropriate 
measure of the escape rate. Secondly, not all escapers from the Hill potential actually 
intersect the obvious surface of section ?/ = 0. Thirdly, the problem is three-dimensional, 
and the visualisation of turnstile dynamics becomes harder; Wiggins' book shows some 
of the complications that arise. 

On the other hand, in the planar Hill problem some results are possible. In particular, 
the analogues of the escape and capture lobes, Li_2(l) and L2,i(l), and their iterates have 
been mapped out at one or two values of the energy (Roy 2000, Simo & Stuchi 2000). 
For small numbers of iterations one obtains fairly simple ovals on the surface of section. 
These are the intersections of the surface of section with the stable and unstable invariant 
manifolds of the Liapounov orbits, i.e. structures like the tube in Figure ^. For higher 
numbers of iterations their structure becomes highly convoluted, and further complicated 
by the fact that, at some intersections, only part of the tube actually intersects the surface. 

Another factor which turnstile dynamics clarifies is the relationship between escape, 
which is our interest here, and temporary capture, which has motivated other studies (e.g. 
Murison 1989). 
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